//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

// path of an electron aimed at the funny cusps/line cusps with gnuplot

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include "ephi.hpp"

#include <vector>

using namespace std;

#define RADIUS 0.15
#define SPACING 0.01
#define WR 0.035
#define CURRENT 200000
#define CHARGEPM 3e-7

int main (int argc, char *argv[])
{
  Statics statics;
  make_polywell_cube (statics, RADIUS, WR, SPACING, CURRENT, 0);
  ElectroDynamics dyn(statics);
  dyn.setDisableFields (true);
  dyn.setDTM (1, 10);

  FILE *pf = fopen ("test24.plot", "w");
  fprintf (pf, "set style data lines\n");
  fprintf (pf, "unset title\n");
  fprintf (pf, "set ticslevel 0\n");
  fprintf (pf, "set xrange [-0.3:0.3]\n");
  fprintf (pf, "set yrange [-0.3:0.3]\n");
  fprintf (pf, "set zrange [-0.3:0.3]\n");
  fprintf (pf, "set view 75, 30;\n");
  fprintf (pf, "set view 0, 0;\n");
  fprintf (pf, "set terminal png nocrop enhanced size 1200,1200\n");
  fprintf (pf, "set output 'test24.png'\n");
  fprintf (pf, "set multiplot\n");
  std::string str = statics.get_gnuplot_script().c_str();
  fprintf (pf, "%s", str.c_str());

#define SPEED 10000

  prec_t ticks = dyn.ticks_for_distance (ElectroDynamics::eV2ms (SPEED, ELECTRON_MASS), 0.5);

  prec_t phi = PREC_PI / 2880;
  prec_t angle = PREC_PI / 4;
  vect3d maindir(cos(angle), sin(angle), 0);
  for ( size_t i = 0; i < 40; i++ )
    {
      vect3d dir(cos(angle), sin(angle), 0);
      dyn.inject_electron (vect3d(0, 0, 0), dir.normal(), SPEED);
      fprintf (pf, "splot '-' notitle\n");
      for ( size_t i = 0; i < 100; i++ )
        {
          for ( size_t j = 0; j < ticks / 100; j++ )
              dyn.step();
          vect3d pos = dyn.get_pos(0);
          fprintf (pf, "%.3e %.3e %.3e\n", prec2double(pos.x), prec2double(pos.y), prec2double(pos.z));
          if ( pos.length() > 0.3 )
              break;
        }
      printf ("Angle is: %g\n", prec2double(360 * (angle - PREC_PI / 4) / (2 * PREC_PI)));
      if ( dyn.get_pos(0).length() >= 0.3 )
          printf ("pass!\n");
      else
          printf ("fail!\n");
      fprintf (pf, "e\n");
      angle += phi;
      dyn.erase (0);
    }

  fprintf (pf, "unset multiplot\n");

  fclose (pf);

  return 0;
}

